#----------------- #
# Plots and tables #
#----------------- #
rm(list = ls())

## Function for later
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}


## data sets
adjacency_matricies <- readRDS("intl_order/adjacency_matricies_replication.rds")
community.blondel <- readRDS("intl_order/international_order_replication.rds")
plot_assortativity <- readRDS("data/plot_assortativity.rds")
load("model_output/permutation_estimates_emulate.rda")
load("model_output/permutate_strength_estimates.rda")
polity.diffusion <- readRDS("data/polity_diffusion_replication.rds")
war.data <- readRDS("data/war_data_replication.rds")
load("model_output/conflict_model_fits.Rda")
pred_values <- readRDS("model_output/pred_values_outgroup.rds")

gc()
## Conflict models 
load("model_output/multi_factor_summary.RData")
load("model_output/full_lpm.Rda")
load("model_output/full_regression.Rda")
load("model_output/btergm_output.Rdata")

## Diffusion models 
load("model_output/re.exec.const.Rda")
load("model_output/re.autoc.Rda")
load("model_output/re.democ.Rda")
load("model_output/re.rest.part.Rda")
load("model_output/re.free.media.Rda")

load("model_output/point.est.re.exec.const.Rda")
load("model_output/point.est.re.autoc.Rda")
load("model_output/point.est.re.democ.Rda")
load("model_output/point.est.re.rest.part.Rda")
load("model_output/point.est.re.free.media.Rda")
intl_groups <- readRDS("intl_order/international_order_replication.rds")
coop_regime <- dplyr::select(
  intl_groups, 
  c(
    ccode, year, interlayer
  )
)

## Assortativity models 
load("model_output/granger_assortatative_regime.RData")
load("model_output/var_assortativity.RData")

polity.diffusion <- mutate(
  polity.diffusion,
  democ = ifelse(polity >= 6, 1, 0),
  autoc = ifelse(polity <= -6, 1, 0),
  v2csgender_ord = as.numeric(as.character(v2csgender_ord)),
  v2csrlgrep_ord = as.numeric(as.character(v2csrlgrep_ord)),
  parcomp = as.numeric(as.character(parcomp)),
  parreg = as.numeric(as.character(parreg)),
  xconst = as.numeric(as.character(xconst)),
  v2mecenefm_ord = as.numeric(as.character(v2mecenefm_ord)),
  women_suffrage = ifelse(v2csgender_ord >= 3, 1, 0),
  religious_freedom = ifelse(v2csrlgrep_ord >= 3, 1, 0),
  free_media = ifelse(v2mecenefm_ord >= 3, 1, 0),
  wdi_gdpcapgr = as.numeric(scale(wdi_gdpcapgr)),
  comptetive_political_participation = ifelse(parcomp == 5, 1, 0),
  unilateral_executive = ifelse(xconst  == 1, 1, 0),
  restricted_political_participation = ifelse(parreg == 4, 1, 0),
  executive_constrainted = ifelse(xconst == 7, 1, 0)
)

## Plot degree
temp_net <- graph_from_adjacency_matrix(adjacency_matricies[[1]], mode = "undirected", weighted = T)
degree_data <- data.frame(
  names = names(strength(temp_net, mode = "total")),
  weighted = strength(temp_net, mode = "total"),
  year = 1816,
  variable = "Degree centrality"
)
eigenvector_data <- data.frame(
  names = names(eigen_centrality(temp_net, directed = F, weights = E(temp_net)$weight, scale = F)$vector),
  weighted = eigen_centrality(temp_net, directed = F, weights = E(temp_net)$weight, scale = F)$vector,
  year = 1816,
  variable = "Eigenvector centrality"
)

for (i in 2:length(adjacency_matricies))
{
  g <- graph_from_adjacency_matrix(adjacency_matricies[[i]], mode = "undirected", weighted = T)
  degree_data <- bind_rows(
    degree_data,
    data.frame(
      names = names(strength(g)),
      weighted = strength(g, mode = "total"),
      year = 1815 + i,
      variable = "Degree centrality"
    )
  )
  
  eigenvector_data <- bind_rows(
    eigenvector_data,
    data.frame(
      names = names(eigen_centrality(g, directed = F, weights = E(g)$weight, scale = F)$vector),
      weighted = eigen_centrality(g, directed = F, weights = E(g)$weight, scale = F)$vector,
      year = 1815 + i,
      variable = "Eigenvector centrality"
    )
  )
}

centrality_appendix_table <- list()
for (i in 1816:2014)
{
  highest_deg <- filter(
    degree_data,
    year == i
  ) %>% arrange(
    -weighted
  ) 
  highest_deg <- highest_deg[1:10,]
  highest_deg$rank_deg = 1:10
  highest_deg <- dplyr::rename(
    highest_deg, 
    state_deg = names, 
    degree = weighted
  ) %>% dplyr::select(
    -c(variable, year)
  )
  highest_eig <- filter(
    eigenvector_data,
    year == i
  ) %>% arrange(
    -weighted
  )
  highest_eig <- highest_eig[1:10,]
  highest_eig$rank_eig = 1:10
  highest_eig <- dplyr::rename(
    highest_eig, 
    state_eig = names, 
    eigenvector = weighted
  ) %>% dplyr::select(
    -c(variable, year)
  )
  output_table_centrality <- bind_cols(
    data.frame(
      year = i
    ),
    highest_deg,
    highest_eig
  )
  rownames(output_table_centrality) <- NULL
  centrality_appendix_table[[length(centrality_appendix_table) + 1]] <- output_table_centrality
}

all_centrality <- do.call(bind_rows, centrality_appendix_table)
all_centrality$eigenvector <- round(all_centrality$eigenvector, 2)
all_centrality <- dplyr::select(
  all_centrality, 
  c(
    year, rank_deg, state_deg, degree, state_eig, eigenvector
  )
)
for (i in 1:ncol(all_centrality))
{
  if (!(i %in% c(1, ncol(all_centrality))))
  {
    if (is.numeric(all_centrality[,i]))
    {
      all_centrality[,i] <- paste0("$", all_centrality[,i], "$&")   
    } else {
      all_centrality[,i] <- paste0(" ", all_centrality[,i], "&")   
    } 
  } else if (i == 1)
  {
    if (is.numeric(all_centrality[,i]))
    {
      all_centrality[,i] <- paste0("$", all_centrality[,i], "$&")   
    } else {
      all_centrality[,i] <- paste0("", all_centrality[,i], "&")   
    }
  } else
  {
    if (is.numeric(all_centrality[,i]))
    {
      all_centrality[,i] <- paste0("$", all_centrality[,i], "$\\\\")   
    } else {
      all_centrality[,i] <- paste0("", all_centrality[,i], "\\\\")   
    }
  }
  
}

all_centrality$eigenvector[seq(10, nrow(all_centrality), 10)] <- paste0(all_centrality$eigenvector[seq(10, nrow(all_centrality), 10)], "\\hline")

write.table(all_centrality, file = "data_visualization/si_table_4.txt", row.names = F, sep = " ")

degree_data$weighted <- range01(degree_data$weighted)
eigenvector_data$weighted <- range01(eigenvector_data$weighted)

degree_data <- bind_rows(
  degree_data
)

powerful_states <- filter(
  degree_data,
  names %in% c("USA", "CHN", "UKG", "FRN", "RUS")
)
some_states <- ggplot(subset(degree_data, !(names %in% powerful_states$names)), aes(x = year, y = weighted, group = names)) + 
  geom_line(alpha = I(1/8)) + 
  facet_wrap(~variable) + 
  theme_bw() + 
  labs(x = "Year", y = "Weighted centrality") + 
  theme(strip.text.x = element_text(size = 18, color = "black"), 
        axis.text = element_text(size = 18,  color = "black"), 
        axis.title = element_text(size = 18,  color = "black"),
        axis.ticks.x = element_blank())

some_states + geom_line(data = subset(powerful_states), aes(x = year, y = weighted, group = names, col = names, linetype = names), linewidth = 1) + 
  scale_colour_brewer(palette = "Dark2") + 
  theme(legend.position = "bottom", legend.title = element_blank(),
        strip.text.x = element_text(size = 18, color = "black"), 
        axis.text = element_text(size = 18,  color = "black"), 
        axis.title = element_text(size = 18,  color = "black"),
        legend.text = element_text(size = 18,  color = "black")) +
  guides(color=guide_legend(nrow = 1, byrow = TRUE))


ggsave("data_visualization/centrality_over_time.png", width = 7, height = 5, units = "in")  

## Permutation data
check_order <- filter(
  community.blondel,
  year == 1980
)

gc()
permutation_estimates$perm_democ         <- permutation_estimates$perm_democ - fixef(point.est.re.democ)[2]
permutation_estimates$perm_autoc         <- permutation_estimates$perm_autoc - fixef(point.est.re.autoc)[2]
permutation_estimates$perm_pol_rest_part <- permutation_estimates$perm_pol_rest_part - fixef(point.est.re.exec.const)[2]
permutation_estimates$perm_free_media    <- permutation_estimates$perm_free_media - fixef(point.est.re.free.media)[2]
permutation_estimates$perm_exec_const    <- permutation_estimates$perm_exec_const - fixef(point.est.re.exec.const)[2]

permutation_strength_estimates$perm_democ         <- permutation_strength_estimates$perm_democ - fixef(point.est.re.democ)[2]
permutation_strength_estimates$perm_autoc         <- permutation_strength_estimates$perm_autoc - fixef(point.est.re.autoc)[2]
permutation_strength_estimates$perm_pol_rest_part <- permutation_strength_estimates$perm_pol_rest_part - fixef(point.est.re.exec.const)[2]
permutation_strength_estimates$perm_free_media    <- permutation_strength_estimates$perm_free_media - fixef(point.est.re.free.media)[2]
permutation_strength_estimates$perm_exec_const    <- permutation_strength_estimates$perm_exec_const - fixef(point.est.re.exec.const)[2]

plot_data <- reshape2::melt(permutation_estimates)
plot_data$type = "Permutation of frequent interactions"
plot_data_strength <- reshape2::melt(permutation_strength_estimates)
apply(permutation_strength_estimates, 2, function(x){sum(x > 0)})/nrow(permutation_strength_estimates)

plot_data_strength$type = "Permutation of CINC"
plot_data <- bind_rows(
  plot_data, 
  plot_data_strength
)
plot_data$variable <- as.character(plot_data$variable)
plot_data$variable[plot_data$variable == "perm_democ"] <- "Democractic"
plot_data$variable[plot_data$variable == "perm_autoc"] <- "Autocractic"
plot_data$variable[plot_data$variable == "perm_pol_rest_part"] <- "Restricted political participation"
plot_data$variable[plot_data$variable == "perm_exec_const"] <- "Constrained executive"
plot_data$variable[plot_data$variable == "perm_free_media"] <- "Free media"

plot_data$variable <- gsub(" ", "\n", plot_data$variable)
plot_data$variable <- factor(
  plot_data$variable, 
  levels = c(
    "Constrained\nexecutive",
    "Free\nmedia",
    "Restricted\npolitical\nparticipation",
    "Restrictions\non\nreligious\npractices",
    "Autocractic",
    "Democractic"
  )
)
plot_data$label <- "Network simulations frequent interactions and political institutions"

ggplot(plot_data, aes(x = variable, y = value, fill = variable)) + 
  geom_jitter(alpha = I(1/20)) + 
  geom_boxplot(size = 1) + 
  theme_bw() + 
  scale_fill_brewer(palette = "Dark2") + 
  theme(legend.position = "none",
        strip.text.x = element_text(size = 20, color = "black"), 
        axis.text = element_text(size = 15,  color = "black"), 
        axis.title = element_text(size = 20,  color = "black")) + 
  geom_hline(yintercept = 0, size = 2, linetype = "dashed") + 
  facet_wrap(~type) + 
  labs(x = "Diffusion over relational ties", y = "Simulated values centered at zero")

ggsave("data_visualization/plot_permutations.png", width = 13, height = 7, units = "in")  

## Modularity over time
ggplot(plot_assortativity, 
       aes(x = year, y = value, col = variable,
           linetype = variable, shape = variable)) + 
  geom_smooth(alpha = I(2/3), span = 0.2, se = F, size = 2.5) + 
  #geom_line(alpha = I(3/10), size = 1.5) +
  theme_bw() + 
  facet_wrap(~title) + 
  scale_color_manual(values = c("#1B9E77", "#E6AB02", "#D95F02")) + 
  labs(x = "", y = "Scaled values") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.background=element_blank(),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 20, color = "black"), 
        axis.text = element_text(size = 20,  color = "black"), 
        axis.title = element_text(size = 20,  color = "black")) + 
  annotate("rect", xmin = 1945, xmax = 1991, ymin = -Inf, ymax = Inf,
           alpha = .2, fill = "#7570B3") + 
  annotate("rect", xmin = 1939, xmax = 1945, ymin = -Inf, ymax = Inf,
           alpha = .5, fill = "#E7298A") +
  annotate("rect", xmin = 1918, xmax = 1939, ymin = -Inf, ymax = Inf,
           alpha = .2, fill = "#7570B3") +
  annotate("rect", xmin = 1914, xmax = 1918, ymin = -Inf, ymax = Inf,
           alpha = .5, fill = "#E7298A") +
  annotate("rect", xmin = 1816, xmax = 1914, ymin = -Inf, ymax = Inf,
           alpha = .2, fill = "#7570B3") +
  annotate("rect", xmin = 1991, xmax = 1992, ymin = -Inf, ymax = Inf,
           alpha = .5, fill = "#E7298A") +
  annotate("rect", xmin = 1992, xmax = 2010, ymin = -Inf, ymax = Inf,
           alpha = .2, fill = "#7570B3") + 
  annotate("text", x =  1864, y = 1.1, label = "Concert of\n Europe", size = 7) + 
  annotate("text", x =  1928.5, y = 1.1, label = "Interwar\n Period", size = 7) + 
  annotate("text", x =  1968, y = 1.1, label = "Cold\n War", size = 7) + 
  annotate("text", x =  2004, y = 1.1, label = "Liberal\n Order", size = 7) + 
  scale_y_continuous(limits = c(0, 1.15), breaks = seq(0, 1, 0.25))

ggsave("data_visualization/modularity_assortativity.png", width = 9, height = 7, units = "in")

## Ego-net data
years <- 1816:2014
index_year <- which(years == 1980)
ego_net_index_year <- adjacency_matricies[index_year]
ego_nets <- lapply(ego_net_index_year, graph_from_adjacency_matrix, mode = "directed", weighted = T)
edgelist <- lapply(ego_nets, as_edgelist, names = TRUE)
edgelist <- lapply(edgelist, data.frame)
edgelist <- do.call(bind_rows, edgelist)
edgelist$weight <- E(ego_nets[[1]])$weight
edgelist <- 
  distinct(
    edgelist,
    X1, X2, .keep_all = T
  )

edgelist_usa <- 
  filter(
    edgelist,
    X1 == "USA" | 
      X2 == "USA"
  ) %>% 
  mutate_all(
    as.character
  )

polity_data <- readxl::read_excel("data/p5v2018d.xls")
polity_data <- filter(
  polity_data, 
  byear <= 1980,
  eyear >= 1980,
  !is.na(polity)
)
polity_data <- unique(polity_data$scode)
edgelist_usa <- unique(c(edgelist_usa$X1, edgelist_usa$X2))
edgelist_usa <- edgelist_usa[edgelist_usa %in% polity_data]
powerful_states <- filter(
  polity.diffusion, 
  year == 1980
) %>% 
  dplyr::select(
    c(
      names, cinc
    )
  )

usa_edgelist <- mutate_all(
  edgelist,
  as.character
) %>% 
  filter(
    X1 %in% edgelist_usa,
    X2 %in% edgelist_usa
  ) %>% 
  as.matrix() 

temp <- usa_edgelist %>% 
  data.frame(.) %>% 
  left_join(
    powerful_states,
    by = 
      c(
        "X1" = "names"
      )
  ) %>% 
  left_join(
    powerful_states,
    by = 
      c(
        "X2" = "names"
      )
  ) %>% 
  mutate(
    check = ifelse(cinc.x > cinc.y, T, F)
  )
#usa_edgelist <- data.frame(usa_edgelist)
#usa_edgelist$weight <- as.numeric(usa_edgelist$weight)


usa <- graph_from_edgelist(usa_edgelist[temp$check == T,1:2], directed = T)
E(usa)$weight <- as.numeric(usa_edgelist[temp$check == T,3])
usa <- igraph::simplify(usa, remove.multiple=T)
colrs <- c("#1B9E77", "#D95F02", "#7570B3")
igraph::vertex_attr(usa, "polity") <- ifelse(V(usa)$name %in% subset(polity.diffusion, year == 1980 & polity >= 6)$names, "Democractic state", ifelse(V(usa)$name %in% subset(polity.diffusion, year == 1980 & polity <= -6)$names, "Autocratic state", "Mixed regime"))

usa_1980 <- filter(
  polity.diffusion,
  year == 1980,
  names == "USA"
)

V(usa)$color <- colrs[as.numeric(as.factor(V(usa)$polity))]
V(usa)$shape <- c("circle", "circle", "circle")[as.numeric(as.factor(V(usa)$polity))]
plot(usa, vertex.label.color = "black", vertex.label.cex=1, vertex.size = 11)
legend(x= -1, 
       y= -0.75, 
       c("Democractic state", "Mixed regime", "Autocratic state"), pch=c(21, 21, 21),
       col = colrs[c(2, 3, 1)], 
       pt.bg = colrs[c(2, 3, 1)], pt.cex=2.5, cex=1.5, bty="n", ncol=1)



temp_1980 <- filter(
  polity.diffusion, 
  year == 1980
) %>% 
  arrange(
    names
  )
plot_community_data <- adjacency_matricies[which(years %in% c(1980))]
plot_community_data <- graph_from_adjacency_matrix(plot_community_data[[1]], mode = "undirected", weighted = T)
plot_community_data <- igraph::simplify(plot_community_data, remove.multiple=T)
igraph::vertex_attr(plot_community_data, "group") <-  temp_1980$membership

#Isolated = which(degree(plot_community_data)==0)
#G2 = igraph::delete.vertices(plot_community_data, Isolated)
V(plot_community_data)$color <- c("#E7298A", "#66A61E", "#E6AB02", "cadetblue", "yellow", "purple")[as.numeric(as.factor(V(plot_community_data)$group))]
plot(plot_community_data, vertex.label.color = "black", vertex.label.cex=1, vertex.size = 11, layout = layout_with_fr)

out_group_pred <- filter(
  pred_values,
  out_group %in% 0:1
) %>% mutate(
  variable = ifelse(out_group == 1, "Out-group relation", "In-group relation"),
  facet_lab = "In- versus out-group relation"
) 

out_group_pred %>% 
  plyr::ddply(
    ~variable,
    summarize,
    est = mean(pred, na.rm = T)
  )

joint_democ_pred <- filter(
  pred_values,
  democ_joint %in% 0:1
) %>% mutate(
  variable = ifelse(democ_joint == 1, "Joint democracies", "Not joint democracies"),
  facet_lab = "Joint versus not joint democracy"
)


joint_democ_pred %>% 
  plyr::ddply(
    ~variable,
    summarize,
    est = mean(pred, na.rm = T)
  )

pred_values_plot <- bind_rows(
  out_group_pred, 
  joint_democ_pred
) %>% mutate(
  variable = factor(
    variable, 
    levels = c(
      "In-group relation",
      "Out-group relation",
      "Joint democracies",
      "Not joint democracies"
    )
  )
)

ggplot(pred_values_plot, aes(x = variable, y = pred, fill = variable)) + 
  geom_jitter(alpha = I(1/15), size = 1.5) + 
  geom_boxplot(size = 1.25) + 
  theme_bw() + 
  scale_colour_brewer(palette = "Dark2") + 
  scale_fill_brewer(palette = "Dark2") + 
  theme(legend.position = "bottom",
        legend.title = element_blank()) + 
  labs(x = "", y = "Bootstrapped predicted value (USA-JPN in 1941)") + 
  facet_wrap(~facet_lab, scales = "free_x") + 
  theme(legend.position = "none",
        legend.text = element_text(size = 15, color = "black"),
        strip.text.x = element_text(size = 15, color = "black"), 
        axis.text = element_text(size = 15,  color = "black"), 
        axis.title = element_text(size = 15,  color = "black")) 

ggsave("data_visualization/democ_out_group.png",  width = 12, height = 7, units = "in")  

#define object to plot
ggplot(subset(all_fits, model != "TERGM"), aes(m = pred, d = result, col = model)) + 
  geom_roc(n.cuts = 50, labels = F, size = 1) + 
  facet_wrap(~dv) + 
  theme_bw() + 
  scale_colour_brewer(palette = "Dark2") + 
  theme(legend.position = "bottom",
        legend.title = element_blank()) + 
  labs(x = "False positive fraction",
       y = "True positive fraction",
       legend.text = element_text(size = 20),
       strip.text.x = element_text(size = 20, color = "black"), 
       axis.text = element_text(size = 20,  color = "black"), 
       axis.title = element_text(size = 20,  color = "black"))


models <- unique(all_fits$model)
all_fits <- subset(all_fits, !is.na(pred))
output_pres <- data.frame()
for (i in 1:length(models))
{
  temp_result <- evalmod(scores = subset(all_fits, model == models[i])$pred, 
                         labels = subset(all_fits, model == models[i])$result)
  output_pres <- data.frame(
    curve = c(rep("ROC", length(temp_result$rocs[[1]]$x)), rep("Precision-recall", length(temp_result$prcs[[1]]$x))),
    x = c(temp_result$rocs[[1]]$x, temp_result$prcs[[1]]$x),
    y = c(temp_result$rocs[[1]]$y, temp_result$prcs[[1]]$y),
    model = models[i],
    label = "Fatal MIDs model fit"
  ) %>% bind_rows(
    output_pres
  )
}
output_pres <- transform(
  output_pres, 
  model = ifelse(model == "Multiplicative\nlatent factor model", "Latent\nfactor model", model),
  curve = factor(curve, levels = c("ROC", "Precision-recall"))
)

roc_plot <- ggplot(subset(output_pres, curve == "ROC"), aes(x = x, y = y, col = model)) + 
  facet_wrap(~curve) + 
  geom_line(size = 1) + 
  scale_colour_brewer(palette = "Dark2") + 
  theme_bw() + 
  theme(legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size = 10),
        strip.text.x = element_text(size = 15, color = "black"), 
        axis.text = element_text(size = 15,  color = "black"), 
        axis.title = element_text(size = 15,  color = "black")) + 
  labs(x = "1 - Specificity", y = "Sensitivity")


proc_plot <- ggplot(subset(output_pres, curve != "ROC"), aes(x = x, y = y, col = model)) + 
  facet_wrap(~curve, scales = "free") + 
  geom_line(size = 1) + 
  scale_colour_brewer(palette = "Dark2") + 
  theme_bw() + 
  theme(legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, color = "black"), 
        axis.text = element_text(size = 15,  color = "black"), 
        axis.title = element_text(size = 15,  color = "black")) + 
  labs(y = "Precision", x = "Recall")


roc_proc_fit <- cowplot::plot_grid(roc_plot, proc_plot, ncol = 2)
ggsave("data_visualization/conflict_model_fit.png", roc_proc_fit, width = 10, height = 5, units = "in")  

## Semiparametric bootstrap 

concert1 <- subset(
  polity.diffusion, 
  year < 1831
) %>% 
  distinct(
    names, membership
  )

concert2 <- subset(
  polity.diffusion, 
  year >= 1910 & 
    year < 1914
) %>% 
  distinct(
    names, membership
  )

interwar <- subset(
  polity.diffusion, 
  year > 1941 &  
    year < 1945
) %>% 
  distinct(
    names, membership
  )

coldwar <- subset(
  polity.diffusion, 
  year >= 1980 &  
    year < 1985
) %>% 
  distinct(
    names, membership
  )


liberal <- subset(
  polity.diffusion, 
  year >= 1991
) %>% 
  distinct(
    names, membership
  )


range <- 20
year_start <- seq(1816, 2014, range)
year_end   <- unique(c(seq(1816 + range - 1, 2014, range), 2014))

cv_time <- list()
for (i in 1:length(year_start))
{
  suppressMessages({
    suppressWarnings({
      vary_time <- data.frame(matrix(NA, nrow = 5, ncol = 4))
      colnames(vary_time) <- c("coef", "mean", "se", "year")
      temp_data <- as.data.frame(subset(war.data, !(year %in% year_start[i]:year_end[i])))
      temp_data <- dplyr::select(
        temp_data,
          c(
            fatal_mid, out_group, democ_joint,
            autoc_joint, major_power_joint, cinc_ratio,
            namea, nameb, dyad, year
            )
       )
      temp_model <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | namea + nameb + dyad  + year, data = temp_data, family = "binomial")
      vary_time[,2:3] <- summary(temp_model, cluster = c("namea", "nameb", "dyad", "year"))$coeftable
      vary_time$coef <- rownames(summary(temp_model)$coeftable)
      vary_time$year <- paste(year_start[i], year_end[i], sep = "-")
      vary_time$significant <- ifelse(summary(temp_model, cluster = c("namea", "nameb", "dyad", "year"))$coeftable[,4] < 0.05, "Statistically\nsignficant", "Not statistically\nsignificant")
      cv_time[[i]] <- vary_time
      rm(temp_model)
    })
  })
}

cv_time <- do.call(bind_rows, cv_time)
cv_time$coef[cv_time$coef == "out_group"] <- "Out-group\nrelation"
cv_time$coef[cv_time$coef == "democ_joint"] <- "Joint\ndemocracy"
cv_time$coef[cv_time$coef == "autoc_joint"] <- "Joint\nautocracy"
cv_time$coef[cv_time$coef == "major_power_joint"] <- "Joint major\n power"
cv_time$coef[cv_time$coef == "cinc_ratio"] <- "Cinc\nratio"
cv_time$lb <- cv_time$mean - 1.96 * cv_time$se
cv_time$ub <- cv_time$mean + 1.96 * cv_time$se

cv_time$year <- factor(cv_time$year, 
                       levels = unique(cv_time$year))
ggplot(cv_time, aes(x = mean, y = reorder(coef, -lb), xmin = lb, xmax = ub, col = significant)) + 
  geom_point() + 
  geom_errorbar(size = 2, width = 0.5) + 
  facet_wrap(~year, ncol = 2) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  theme(legend.position = "bottom", 
        legend.title = element_blank(),
        strip.text.x = element_text(size = 20, color = "black"),
        axis.text = element_text(size = 20,  color = "black"),
        axis.title = element_text(size = 20,  color = "black"),
        legend.text = element_text(size = 20)) + 
  labs(x = "Cross validated estimates\n(leaving out time)", 
       y = "Dropped time interval")

ggsave("data_visualization/cv_data.png", width = 13, height = 18, units = "in")  


layers_to_check <- data.frame(
  expand.grid(
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1),
    c(0, 1)
  )
) 


layers_to_check <- layers_to_check[rowSums(layers_to_check) > 0,]
colnames(layers_to_check) <- c("atop_defense",
                               "atop_offense",
                               "atop_consul",
                               "atop_neutral",
                               "bla", 
                               "dipl", 
                               "igo",
                               "pta", 
                               "dcaV2", 
                               "dcaV1"
)
layers_to_check[is.na(layers_to_check)] <- 0

layers_to_check <- filter(
  layers_to_check,
  !(dcaV1 == 1 & dcaV2 == 1),
  atop_defense + atop_offense + igo > 0
)

layers_to_check$layer_combo <- 1:nrow(layers_to_check)



formula_list <- data.frame(
  ## Fixed 
  formula_downweight = 1:12,
  formula = c(
    "tie ~\ncontig + year",
    "tie ~\ncontig + year\n+ major_power",
    "tie ~\ncontig + year +\nyear_sqr + major_power",
    "tie ~\ncontig  | year",
    "tie ~\ncontig + joint_region\n| year",
    "tie ~\ncontig\n|  year + cowc1 + cowc2",
    "tie ~\ncontig + joint_region\n| year + cowc1 + cowc2",
    "tie ~\ncontig + major_power\n| year + cowc1 + cowc2",
    "tie ~\ncontig + joint_region + major_power\n| year + cowc1 + cowc2",
    "tie ~\ncontig\n| year + cowc1 + cowc2 +\nregion1 + region2",
    "tie ~\ncontig + major_power\n| year + cowc1 +\ncowc2 + region1 + region2",
    "No\nweight")
)

outcome_data_1  <- readRDS("loop_mods/outcome_data_1.rds")
outcome_data_2  <- readRDS("loop_mods/outcome_data_2.rds")
outcome_data_3  <- readRDS("loop_mods/outcome_data_3.rds")
outcome_data_4  <- readRDS("loop_mods/outcome_data_4.rds")
outcome_data_5  <- readRDS("loop_mods/outcome_data_5.rds")
outcome_data_6  <- readRDS("loop_mods/outcome_data_6.rds")
outcome_data_7  <- readRDS("loop_mods/outcome_data_7.rds")
outcome_data_8  <- readRDS("loop_mods/outcome_data_8.rds")
outcome_data_9  <- readRDS("loop_mods/outcome_data_9.rds")
outcome_data_10 <- readRDS("loop_mods/outcome_data_10.rds")
outcome_data_11 <- readRDS("loop_mods/outcome_data_11.rds")
outcome_data_12 <- readRDS("loop_mods/outcome_data_12.rds")

avg_outcome_data_1  <- readRDS("loop_mods/avg_outcome_data_1.rds")
avg_outcome_data_2  <- readRDS("loop_mods/avg_outcome_data_2.rds")
avg_outcome_data_3  <- readRDS("loop_mods/avg_outcome_data_3.rds")
avg_outcome_data_4  <- readRDS("loop_mods/avg_outcome_data_4.rds")
avg_outcome_data_5  <- readRDS("loop_mods/avg_outcome_data_5.rds")
avg_outcome_data_6  <- readRDS("loop_mods/avg_outcome_data_6.rds")
avg_outcome_data_7  <- readRDS("loop_mods/avg_outcome_data_7.rds")
avg_outcome_data_8  <- readRDS("loop_mods/avg_outcome_data_8.rds")
avg_outcome_data_9  <- readRDS("loop_mods/avg_outcome_data_9.rds")
avg_outcome_data_10 <- readRDS("loop_mods/avg_outcome_data_10.rds")
avg_outcome_data_11 <- readRDS("loop_mods/avg_outcome_data_11.rds")
avg_outcome_data_12 <- readRDS("loop_mods/avg_outcome_data_12.rds")

outcome_data <- bind_rows(
  outcome_data_1,
  outcome_data_2,
  outcome_data_3,
  outcome_data_4,
  outcome_data_5,
  outcome_data_6,
  outcome_data_7,
  outcome_data_8,
  outcome_data_9,
  outcome_data_10,
  outcome_data_11,
  outcome_data_12
) |> distinct(
  layer_combo, formula_downweight, .keep_all = T
) |> transform(
  weight = "Manhattan distance",
  weight_num = 1
)

avg_outcome_data <- bind_rows(
  avg_outcome_data_1,
  avg_outcome_data_2,
  avg_outcome_data_3,
  avg_outcome_data_4,
  avg_outcome_data_5,
  avg_outcome_data_6,
  avg_outcome_data_7,
  avg_outcome_data_8,
  avg_outcome_data_9,
  avg_outcome_data_10,
  avg_outcome_data_11,
  avg_outcome_data_12
) |> distinct(
  layer_combo, formula_downweight, .keep_all = T
) |> transform(
  weight = "Average distance",
  weight_num = 0
) |> dplyr::rename(
  out_group = out_group_fmids,
  p_value = out_group_mids
)

cor(avg_outcome_data$out_group, outcome_data$out_group_fmids)

colnames(outcome_data)[1:2] <- c("out_group", "p_value")
outcome_data <- bind_rows(
  outcome_data,
  avg_outcome_data
)


outcome_data <- left_join(
  outcome_data, 
  layers_to_check
) |> left_join(
  formula_list
)

outcome_vars    <- lm(out_group ~ atop_defense + atop_offense + atop_neutral + atop_consul + bla + dipl + igo + pta + dcaV1 + dcaV2, data = outcome_data)
outcome_p_value <- lm(p_value ~ atop_defense + atop_offense + atop_neutral + atop_consul + bla + dipl + igo + pta + dcaV1 + dcaV2, data = outcome_data)

form_outcome_vars    <- lm(out_group ~ formula, data = outcome_data)
form_outcome_p_value <- lm(p_value ~ formula, data = outcome_data)


weight_outcome_vars    <- lm(out_group ~ weight_num, data = outcome_data)
weight_outcome_p_value <- lm(p_value ~ weight_num, data = outcome_data)


plot_results <- data.frame(
  mean  = c(
    coef(outcome_vars)[-1] + coef(outcome_vars)[1],
    coef(outcome_p_value)[-1] + coef(outcome_p_value)[1],
    coef(form_outcome_vars)[1],
    coef(form_outcome_vars)[-1] + coef(form_outcome_vars)[1],
    coef(form_outcome_p_value)[1],
    coef(form_outcome_p_value)[-1] + coef(form_outcome_p_value)[1],
    coef(weight_outcome_vars)[1],
    coef(weight_outcome_vars)[-1] + coef(weight_outcome_vars)[1],
    coef(weight_outcome_p_value)[1],
    coef(weight_outcome_p_value)[-1]  + coef(weight_outcome_p_value)[1]
  ),
  names =  c(
    names(coef(outcome_vars))[-1],
    names(coef(outcome_p_value))[-1],
    names(coef(form_outcome_vars)),
    names(coef(form_outcome_p_value)),
    "Average\ndistance",
    names(coef(weight_outcome_vars))[2],
    "Average\ndistance",
    names(coef(weight_outcome_p_value))[2]
  ),
  se = c(
    summary(outcome_vars)$coefficient[-1,2],
    summary(outcome_p_value)$coefficient[-1,2],
    summary(form_outcome_vars)$coefficient[,2],
    summary(form_outcome_p_value)$coefficient[,2],
    summary(weight_outcome_vars)$coefficient[,2],
    summary(weight_outcome_p_value)$coefficient[,2]
  ),
  outcome = c(rep("Coefficient effect", 10), rep("p-value", 10),
              rep("Coefficient effect", 12), rep("p-value", 12),
              rep("Coefficient effect", 2), rep("p-value", 2)),
  model = c(rep("Layer", 20), rep("Formula", 24), rep("Weight", 4))
)
plot_results$lb <- plot_results$mean + -1.96 * plot_results$se
plot_results$ub <- plot_results$mean + 1.96 * plot_results$se

plot_results <- transform(
  plot_results,
  names = gsub("formula", "", names),
  label = paste0(model, " ", outcome)
)


plot_results <- transform(
  plot_results,
  names = factor(
    names,
    levels = c(
      "atop_defense",
      "atop_offense",
      "atop_consul",
      "atop_neutral",
      "bla", 
      "dipl", 
      "igo",
      "pta", 
      "dcaV2", 
      "dcaV1",
      "(Intercept)",
      unique(outcome_data$formula)[order(nchar(unique(unique(outcome_data$formula))))][-1],
      "Average\ndistance",
      "weight_num"
    ),
    labels = c(
      "Defensive\nalliance", "Offensive\nalliance", "Consulation\nalliance", 
      "Neutral\nalliance",
      "BLA", "Ambass-\nador", "IGO", 
      "PTA", "DCA (V1)", "DCA (V2)", "No weight",
      unique(outcome_data$formula)[order(nchar(unique(unique(outcome_data$formula))))][-1],
      "Average\ndistance",
      "Manhattan\ndistance"
    ))
)


layers_plot_eff <- ggplot(plot_results, aes(x = mean, xmin = lb, xmax = ub, y = reorder(names, mean), col = model)) + 
  geom_errorbar(width = 0.25, size = 1) + 
  facet_grid(model~outcome, scales = "free", space = "free_y") + 
  theme_bw() + 
  scale_color_manual(values = c("#7570B3", "#E7298A", "#66A61E")) + 
  theme(legend.position = "none",
        axis.text = element_text(size = 12,  color = "black"),
        legend.title = element_blank(),
        legend.text = element_text(size = 10),
        strip.text = element_text(size = 12, color = "black"), 
        axis.title = element_text(size = 12,  color = "black")) + 
  labs(x = "", y = "")


ggsave("data_visualization/layer_effect.png", layers_plot_eff, width = 12, height = 15, units = "in")  

layer_combo_plot <- reshape2::melt(
  layers_to_check,
  id = "layer_combo"
) |> transform(
  variable = factor(
    variable,
    levels = c(
      "atop_defense",
      "atop_offense",
      "atop_consul",
      "atop_neutral",
      "bla", 
      "dipl", 
      "igo",
      "pta", 
      "dcaV2", 
      "dcaV1"
    ),
    labels = c(
      "Defensive\nalliance", "Offensive\nalliance", "Consulation\nalliance", 
      "Nuetral\nalliance",
      "BLA", "Ambassador", "IGO", 
      "PTA", "DCA (V1)", "DCA (V2)"
    )),
  layer_combo = factor(layer_combo, levels = 1:6702),
  value = ifelse(value == 1, "Layer\nincluded", "Layer\nnot included"),
  label = "Alternative input layers"
)

ggplot(layer_combo_plot, aes(x = layer_combo, y = variable, fill = value)) + 
  geom_tile(colour = "white") + 
  scale_fill_brewer(palette = "Set1") + 
  labs(x = "Layer combination", y = "") +
  theme_bw() + 
  facet_wrap(~label) + 
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6,  color = "black"),
        legend.title = element_blank(),
        legend.text = element_text(size = 10),
        strip.text.x = element_text(size = 12, color = "black"), 
        axis.text.y = element_text(size = 10,  color = "black"), 
        axis.title = element_text(size = 12,  color = "black")) + 
  scale_x_discrete(breaks = function(x){x[c(F, F, F, F, F, F, F, F, F, T)]})

ggsave("data_visualization/layer_combo.png", width = 7, height = 5, units = "in")  

outcome_data$label <- "Alternative conflict\nestimates by weight formula"
outcome_data$forms <- factor(outcome_data$formula, levels = unique(outcome_data$formula)[order(nchar(unique(unique(outcome_data$formula))))])
coefs_alt_weight <- ggplot(outcome_data, aes(y = formula, x = out_group, fill = formula)) + 
  geom_boxplot() + 
  scale_fill_brewer(palette = "Set3") + 
  geom_vline(xintercept = 0, col = "black", linetype = "dashed", size = 1) + 
  theme_bw() + 
  facet_grid(weight~label) + 
  labs(x = "Coefficient estimates", y = "") + 
  theme(legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size = 20),
        strip.text = element_text(size = 18, color = "black"), 
        axis.text.y = element_text(size = 9,  color = "black"), 
        axis.text.x = element_text(size = 20,  color = "black"), 
        axis.title = element_text(size = 20,  color = "black"))


outcome_data <- transform(
  outcome_data,
  sig_pos = ifelse(out_group > 0 & p_value <= 0.1, "Statistically\nsignficant", "Not statistically\nsignficant")
)


alt_layers_pvales <- plyr::ddply(
  outcome_data, 
  .(formula, weight), 
  summarize,
  sig_pos = sum(out_group > 0 & p_value <= 0.1),
  not_pos = sum(out_group < 0) + sum(p_value > 0.1)
) |> transform(
  total = sig_pos + not_pos
) |> transform(
  sig_pos = sig_pos/total,
  not_pos = not_pos/total
) |> dplyr::select(
  -total
)|> reshape2::melt(
  id = c("formula", "weight")
) |> transform(
  variable = ifelse(variable == "sig_pos", "Statistically\nsignificant", "Not statistically\nsignificant"),
  formula = factor(formula, levels = unique(formula)[order(nchar(unique(unique(formula))))])
)

plyr::ddply(
  outcome_data, 
  .(formula, weight), 
  summarize,
  median_value = median(out_group),
  lq = quantile(out_group, 0.25),
  uq = quantile(out_group, 0.75)
) %>% arrange(median_value)

plyr::ddply(
  outcome_data, 
  .(formula, weight), 
  summarize,
  pos_vals = sum(out_group > 0)/(nrow(avg_outcome_data_1))
) %>% arrange(pos_vals)

plyr::ddply(
  outcome_data, 
  .(formula, weight), 
  summarize,
  pos_vals = sum(out_group > 0 & p_value <= 0.1)/(nrow(avg_outcome_data_1))
) %>% arrange(pos_vals)

alt_layers_pvales$wrap_lab <- "p-values by\nalternative weight formula"
alt_pvalues <- ggplot(alt_layers_pvales, aes(x = value, y = formula, fill = variable)) + 
  geom_col() + 
  scale_fill_brewer(palette = "Set1") + 
  theme_bw() + 
  facet_grid(weight~wrap_lab) + 
  labs(x = "p-values < 0.05 (one-sided)", y = "") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 18, color = "black"), 
        axis.text = element_text(size = 20,  color = "black"), 
        axis.title = element_text(size = 20,  color = "black")) + 
  scale_y_discrete(labels = NULL)

plot_alts <- cowplot::plot_grid(coefs_alt_weight, alt_pvalues, ncol = 2, align = "h", axis = "bt", rel_widths = c(1.4, 1))
ggsave("data_visualization/alt_coef_estimates.png", plot_alts, width = 13, height = 14, units = "in")  

outcome_data %>% plyr::ddply(
  ~forms,
  summarize,
  median_val = sum(out_group > 0)/(nrow(avg_outcome_data_1) * 2)
) %>% arrange(
  median_val
)
## Create tables 
set.seed(1)
war.data$fatal_mids_num <- war.data$fatal_mid
war.data$fatal_mids_num <- as.numeric(as.character(war.data$fatal_mids_num))
war.data$year_num <- war.data$year - mean(war.data$year)
war.data$state_count <- as.numeric(war.data$state_count)


training_rows <- sample(nrow(war.data), round(0.75 * nrow(war.data)))
war.data.training <- war.data[training_rows, ]
war.data.testing  <- war.data[-training_rows, ]

extract.brms <- function(model,
                         ...) {
  s <- model
  coefficient.names <- rownames(s$fixed)
  mean <- s$fixed[,1]
  se   <- s$fixed[,2]
  pval <- ifelse(s$fixed[,3] > 0 | s$fixed[,4] < 0, 0.049, 0.51)
  
  coefficient.names <- gsub("1", "", coefficient.names)
  coefficient.names <- gsub("\\[\\[i\\]\\]", "", coefficient.names)
  keep_vars <- which(grepl("memory|poly", coefficient.names) == F)
  coefficient.names[coefficient.names == "Intercept"] <- "Constant"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  coefficient.names[coefficient.names == "unaffiliated_joint"] <- "Joint unaffiliated"
  coefficient.names[coefficient.names == "lagged_fatal_mid"] <- "lagged fatal MID"
  coefficient.names[coefficient.names == "out_group"] <- "Out-group relation"
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = mean[keep_vars],
    se = se[keep_vars],
    pvalues = pval[keep_vars]
  )
  return(tr)
}


extract.btergm <- function(model,
                           ...) {
  s <- summary(model)
  coefficient.names <- rownames(s)
  mean <- s[,1]
  lb <- s[,3]
  ub <- s[,4]
  coefficient.names <- gsub("edgecov\\.", "", coefficient.names)
  coefficient.names <- gsub("\\[\\[i\\]\\]", "", coefficient.names)
  keep_vars <- which(grepl("memory|time", coefficient.names) == F)
  coefficient.names[coefficient.names == "edges"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-group relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  coefficient.names[coefficient.names == "joint_autoc_sensitivity"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "nodecov.state_count"] <- "Count of states"
  #coefficient.names[coefficient.names == "gwdeg.fixed.1"] <- "Geometrically weighted degree"
  #coefficient.names[coefficient.names == "gwdeg.fixed.1"] <- "Geometrically weighted degree"
  coefficient.names[coefficient.names == "joint_democ_sensitivity"] <- "Joint democracy"
  tr <- createTexreg(
    coef.names = coefficient.names[keep_vars],
    coef = mean[keep_vars],
    ci.low = lb[keep_vars],
    ci.up = ub[keep_vars]
  )
  return(tr)
}

extract.fixest <- function(model,
                           include.deviance = TRUE,
                           include.nobs = TRUE,
                           ...) {
  require(fixest)
  s <- summary(model, cluster = c("year", "dyad", "ccode1", "ccode2"))
  coefficient.names <- rownames(s$coeftable)
  coefficient.names[coefficient.names == "edges"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-group relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coeftable[, 1]
  se <- s$coeftable[, 2]
  pval <- s$coeftable[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$nobs
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

extract.glm <- function(model,
                        include.deviance = TRUE,
                        include.nobs = TRUE,
                        ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names[coefficient.names == "IIntercept)"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-group relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$df[2]
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

## hierarchical coefficients
extract.re.model <- function(model,
                             include.loglik = TRUE,
                             include.deviance = TRUE,
                             include.nobs = TRUE,
                             ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names[coefficient.names == "prop_autoc_mimic"] <- "Autocratic exposure"
  coefficient.names[coefficient.names == "prop_democ_mimic"] <- "Democratic exposure"
  coefficient.names[coefficient.names == "prop_exec_constraint"] <- "Executive constraint exposure"
  coefficient.names[coefficient.names == "prop_religious_freedom_mimic"] <- "Religious freedom exposure"
  coefficient.names[coefficient.names == "prop_rest_pol_part"] <- "Restricted political participation exposure"
  coefficient.names[coefficient.names == "prop_gov_censorship_mimic"] <- "Free media exposure"
  
  
  
  coefficient.names[coefficient.names == "autocratic_neighbor"] <- "Autocratic neighbor (%)"
  coefficient.names[coefficient.names == "exec_constraint_neighbor"] <- "Executive constraint neighbor  (%)"
  coefficient.names[coefficient.names == "democratic_neighbor"] <- "Democratic neighbor  (%)"
  coefficient.names[coefficient.names == "religious_freedom_neighbor"] <- "Religious freedom neighbor (%)"
  coefficient.names[coefficient.names == "rest_pol_part_neighbor"] <- "Restricted political participation neighbor (%)"
  coefficient.names[coefficient.names == "gov_censorship_neighbor"] <- "Free media neighbor (%)"
  
  coefficient.names[coefficient.names == "(Intercept)"] <- "Constant"
  
  coefficient.names[coefficient.names == "civil_war"] <- "Civil war"
  coefficient.names[coefficient.names == "gle_cgdpc"] <- "GDPC"
  coefficient.names[coefficient.names == "wdi_gdpcapgr"] <- "Annual growth GDPC"
  coefficient.names[coefficient.names == "unaffiliated"] <- "Isolate"
  
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.loglik == TRUE) {
    lik <- logLik(model)
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- nrow(model@frame)
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg( 
    coef.names = coefficient.names[!grepl("poly", coefficient.names)], 
    coef = co[!grepl("poly", coefficient.names)],
    se = se[!grepl("poly", coefficient.names)],
    pvalues = pval[!grepl("poly", coefficient.names)],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  
  return(tr)
}

suppressMessages({suppressWarnings({texreg_lfm     <- extract.brms(multi_sum)})})
suppressMessages({suppressWarnings({texreg_btergm <- extract.btergm(btergm_output)})})
suppressMessages({suppressWarnings({texreg_lpm <- extract.fixest(full_lpm)})})
suppressMessages({suppressWarnings({texreg_logit <- extract.fixest(full_regression)})})


texreg(list(texreg_logit, texreg_lpm, texreg_btergm, texreg_lfm), stars = 0.05, digits = 3)
table_3 <- texreg(list(texreg_logit, texreg_lpm, texreg_btergm, texreg_lfm), stars = 0.05, digits = 3)
write.table(table_3, file = "data_visualization/table_3.txt")
texreg_autoc                     <- extract.re.model(re.autoc)
texreg_democ                     <- extract.re.model(re.democ)
texreg_exec_const                <- extract.re.model(re.exec.const)
texreg_free_media                <- extract.re.model(re.free.media)
texreg_restricted_political_part <- extract.re.model(re.rest.part)

texreg(list(texreg_exec_const, texreg_restricted_political_part, 
            texreg_free_media,
            texreg_autoc, texreg_democ), stars = 0.05, digits = 2, custom.model.names = 
         c("Constrained executive", "Restricted political participation",
            "Free media",
           "Autocratic", "Democratic"), 
       reorder.coef = c(2, 7, 9, 11, 13, 
                        3, 8, 10, 12, 14, 
                        4, 5, 6, 1))

table_2 <- texreg(list(texreg_exec_const, texreg_restricted_political_part, 
            texreg_free_media,
            texreg_autoc, texreg_democ), stars = 0.05, digits = 2, custom.model.names = 
         c("Constrained executive", "Restricted political participation",
           "Free media",
           "Autocratic", "Democratic"), 
       reorder.coef = c(2, 7, 9, 11, 13, 
                        3, 8, 10, 12, 14, 
                        4, 5, 6, 1))
write.table(table_2, file = "data_visualization/table_2.txt")
## hierarchical coefficients
extract.var <- function(model, link,
                        include.loglik = TRUE,
                        ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$varresult[[link]]$coefficients)
  coefficient.names <- paste0(gsub("ingroup\\.l", "Cooperative community modularity (lag ", coefficient.names))
  coefficient.names <- paste0(gsub("assortatative_regime\\.l", "Same-regime assortativity (lag ", coefficient.names))
  coefficient.names[coefficient.names == "const"] <- "Constant"
  coefficient.names[coefficient.names == "trend"] <- "Trend line"
  coefficient.names[grepl("lag", coefficient.names)] <- paste0(coefficient.names[grepl("lag", coefficient.names)], ")")
  
  
  co <- s$varresult[[link]]$coefficients[, 1]
  se <- s$varresult[[link]]$coefficients[, 2]
  pval <- s$varresult[[link]]$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.loglik == TRUE) {
    lik <- logLik(model)
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  
  
  tr <- createTexreg( 
    coef.names = coefficient.names[!grepl("poly", coefficient.names)], 
    coef = co[!grepl("poly", coefficient.names)],
    se = se[!grepl("poly", coefficient.names)],
    pvalues = pval[!grepl("poly", coefficient.names)],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  
  return(tr)
}

modularity_var <- extract.var(var1, link = 1)
assortativity_var <- extract.var(var1, link = 2)


table_4 <- texreg(list(modularity_var, assortativity_var), stars = 0.05, digits = 2, custom.model.names = 
         c("Modularity of cooperative communities", "Same regime ties")
)

write.table(table_4, file = "data_visualization/table_4.txt", row.names = F, sep = " ")

var_models <- list.files("model_output", pattern = "var_assortativity_lag_")
var_models <- var_models[order(as.numeric(gsub('\\D+','', var_models)))]
var_models_loop <- list()


extract.var <- function(model, link,
                        include.loglik = TRUE,
                        ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$varresult[[link]]$coefficients)
  coefficient.names <- paste0(gsub("ingroup\\.l", "Modular. (lag ", coefficient.names))
  coefficient.names <- paste0(gsub("assortatative_regime\\.l", "Assort. (lag ", coefficient.names))
  coefficient.names[coefficient.names == "const"] <- "Constant"
  coefficient.names[coefficient.names == "trend"] <- "Trend line"
  coefficient.names[grepl("lag", coefficient.names)] <- paste0(coefficient.names[grepl("lag", coefficient.names)], ")")
  
  
  co <- s$varresult[[link]]$coefficients[, 1]
  se <- s$varresult[[link]]$coefficients[, 2]
  pval <- s$varresult[[link]]$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.loglik == TRUE) {
    lik <- logLik(model)
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  
  
  tr <- createTexreg( 
    coef.names = coefficient.names[!grepl("poly", coefficient.names)], 
    coef = co[!grepl("poly", coefficient.names)],
    se = se[!grepl("poly", coefficient.names)],
    pvalues = pval[!grepl("poly", coefficient.names)],
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  
  return(tr)
}

for (i in 1:length(var_models))
{
  load(paste0("model_output/", var_models[i]))
  var_models_loop[[length(var_models_loop) + 1]] <- extract.var(var1, link = 1)
  var_models_loop[[length(var_models_loop) + 1]] <- extract.var(var1, link = 2)
}


si_table_11 <- texreg(
  list(
    var_models_loop[[1]], var_models_loop[[2]],
    var_models_loop[[3]], var_models_loop[[4]],
    var_models_loop[[5]], var_models_loop[[6]],
    var_models_loop[[7]], var_models_loop[[8]],
    #var_models_loop[[9]], var_models_loop[[10]],
    var_models_loop[[11]], var_models_loop[[12]],
    var_models_loop[[13]], var_models_loop[[14]],
    var_models_loop[[15]], var_models_loop[[16]],
    var_models_loop[[17]], var_models_loop[[18]],
    var_models_loop[[19]], var_models_loop[[20]]
       ), 
  stars = 0.05, digits = 2, custom.model.names = 
         rep(c("Mod. ", "Reg."), 9)
)
write.table(si_table_11, file = "data_visualization/si_table_11.txt", row.names = F, sep = " ")


dyad_data <- create_dyadyears() %>% 
  add_nmc() %>% 
  add_peace_years() %>% 
  add_cow_mids() %>% 
  add_cow_majors() %>%
  add_democracy() %>% 
  transform(
    cinc_ratio = if_else(cinc1 > cinc2, (cinc1 - cinc2)/(cinc1 + cinc2), (cinc2 - cinc1)/(cinc1 + cinc2)),
    democ_joint = if_else(polity21 >= 6 & polity22 >= 6, 1, 0),
    autoc_joint = if_else(polity21 <= -6 & polity22 <= -6, 1, 0),
    dyad = paste0(ccode1, "-", ccode2),
    mid = cowmidongoing, 
    fatal_mid = if_else(cowmidongoing == 1 & fatality > 0, 1, 0),
    major_power_joint = if_else(cowmaj1 == 1 & cowmaj2 == 1, 1, 0)
  ) %>% left_join(
    coop_regime, 
    by = c(
      "ccode1" = "ccode",
      "year"
    )
  ) %>% left_join(
    coop_regime, 
    by = c(
      "ccode2" = "ccode",
      "year"
    )
  ) %>% mutate(
    out_group = if_else(interlayer.x != interlayer.y, 1, 0),
    ccode1a = if_else(ccode1 < ccode2, ccode1, ccode2),
    ccode2a = if_else(ccode1 > ccode2, ccode1, ccode2),
    ccode1 = ccode1a,
    ccode2, ccode2a
  ) %>% distinct(
    ccode1, ccode2, year, .keep_all = T
  )

extract.fixest <- function(model,
                           include.deviance = TRUE,
                           include.nobs = TRUE,
                           ...) {
  require(fixest)
  s <- summary(model, cluster = c("ccode1", "ccode2", "dyad", "year"))
  coefficient.names <- rownames(s$coeftable)
  coefficient.names[coefficient.names == "edges"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-group relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coeftable[, 1]
  se <- s$coeftable[, 2]
  pval <- s$coeftable[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$nobs
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

load("model_output/full_regression_mid.rda")
load("model_output/simple_regression_mid.rda")
load("model_output/full_glm.rda")
load("model_output/full_lm.rda")
load("model_output/full_glm_mid.rda")
load("model_output/full_lm_mid.rda")
load("model_output/btergm_output_mids.RData")     
load("model_output/full_regression_no_imputes.rda")
load("model_output/full_lpm_no_imputes.rda")
load("model_output/full_regression_no_imputes_mid.rda")
load("model_output/full_lpm_no_imputes_mid.rda")

texreg_lpm_mid                        <- extract.fixest(full_lpm_mid)
texreg_logit_mid                      <- extract.fixest(full_regression_mid)
texreg_lpm_no_fe_fatal                <- extract.glm(full_glm)
texreg_logit_no_fe_fatal              <- extract.glm(full_lm)
texreg_lpm_no_fe_mid                  <- extract.glm(full_glm_mid)
texreg_logit_no_fe_mid                <- extract.glm(full_lm_mid)
texreg_btergm_mid                     <- extract.btergm(btergm_output_mids)
texreg_full_regression_no_imputes     <- extract.fixest(full_regression_no_imputes)
texreg_full_lpm_no_imputes            <- extract.fixest(full_lpm_no_imputes)
texreg_full_regression_no_imputes_mid <- extract.fixest(full_regression_no_imputes_mid)
texreg_full_lpm_no_imputes_mid        <- extract.fixest(full_lpm_no_imputes_mid)


texreg(list(texreg_lpm_mid, texreg_logit_mid, texreg_lpm_no_fe_mid, texreg_logit_no_fe_mid, 
            texreg_btergm_mid, texreg_full_regression_no_imputes_mid, texreg_full_lpm_no_imputes_mid,
            texreg_lpm_no_fe_fatal, texreg_logit_no_fe_fatal,
            texreg_full_regression_no_imputes, texreg_full_lpm_no_imputes), stars = 0.05, digits = 3)


alt_confl <- texreg(list(
  texreg_lpm_mid, texreg_logit_mid, texreg_lpm_no_fe_mid, texreg_logit_no_fe_mid, 
            texreg_btergm_mid, texreg_full_regression_no_imputes_mid, texreg_full_lpm_no_imputes_mid,
            texreg_lpm_no_fe_fatal, texreg_logit_no_fe_fatal,
            texreg_full_regression_no_imputes, texreg_full_lpm_no_imputes), stars = 0.05, digits = 2)

write.table(alt_confl, file = "data_visualization/si_table_7.txt")

load("model_output/lm.point.est.re.democ.rda")
load("model_output/lm.point.est.re.autoc.rda")
load("model_output/lm.point.est.re.rest.part.rda")
load("model_output/lm.point.est.re.free.media.rda")
load("model_output/lm.point.est.re.exec.const.rda")

load("model_output/biv.democ.rda")
load("model_output/biv.autoc.rda")
load("model_output/biv.rest.part.rda")
load("model_output/biv.free.media.rda")
load("model_output/biv.exec.const.rda")

load("model_output/glm.point.est.re.democ.rda")
load("model_output/glm.point.est.re.autoc.rda")
load("model_output/glm.point.est.re.rest.part.rda")
load("model_output/glm.point.est.re.free.media.rda")
load("model_output/glm.point.est.re.exec.const.rda")

coef_data_alt_spec <- data.frame(
  model = c(
    rep("Biv. logit", 10),
    rep("Linear probability model", 40),
    rep("Logistic regression", 40)
  ),
  dv = c(
    rep("Democratic", 2),
    rep("Autocratic", 2),
    rep("Executive\nconstraints", 2),
    rep("Free\nmedia", 2),
    rep("Restricted\npolical\nparticipation", 2),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8)
  ),
  est = c(
    coef(biv.democ),
    coef(biv.autoc),
    coef(biv.exec.const),
    coef(biv.free.media),
    coef(biv.rest.part),
    
    coef(lm.point.est.re.democ),
    coef(lm.point.est.re.autoc),
    coef(lm.point.est.re.exec.const),
    coef(lm.point.est.re.free.media),
    coef(lm.point.est.re.rest.part),
    
    coef(glm.point.est.re.democ),
    coef(glm.point.est.re.autoc),
    coef(glm.point.est.re.exec.const),
    coef(glm.point.est.re.free.media),
    coef(glm.point.est.re.rest.part)
  ),
  names = c(
   names(coef(biv.democ)),
   names(coef(biv.autoc)),
   names(coef(biv.exec.const)),
   names(coef(biv.free.media)),
   names(coef(biv.rest.part)),
   
   names(coef(lm.point.est.re.democ)),
   names(coef(lm.point.est.re.autoc)),
   names(coef(lm.point.est.re.exec.const)),
   names(coef(lm.point.est.re.free.media)),
   names(coef(lm.point.est.re.rest.part)),
   
   names(coef(glm.point.est.re.democ)),
   names(coef(glm.point.est.re.autoc)),
   names(coef(glm.point.est.re.exec.const)),
   names(coef(glm.point.est.re.free.media)),
   names(coef(glm.point.est.re.rest.part))
  ),
  se = c(
    summary(biv.democ)$coefficients[,2],
    summary(biv.autoc)$coefficients[,2],
    summary(biv.exec.const)$coefficients[,2],
    summary(biv.free.media)$coefficients[,2],
    summary(biv.rest.part)$coefficients[,2],
    
    summary(lm.point.est.re.democ)$coefficients[,2],
    summary(lm.point.est.re.autoc)$coefficients[,2],
    summary(lm.point.est.re.exec.const)$coefficients[,2],
    summary(lm.point.est.re.free.media)$coefficients[,2],
    summary(lm.point.est.re.rest.part)$coefficients[,2],
    
    summary(glm.point.est.re.democ)$coefficients[,2],
    summary(glm.point.est.re.autoc)$coefficients[,2],
    summary(glm.point.est.re.exec.const)$coefficients[,2],
    summary(glm.point.est.re.free.media)$coefficients[,2],
    summary(glm.point.est.re.rest.part)$coefficients[,2]
  )
)

coef_data_alt_spec <- filter(
  coef_data_alt_spec,
  !grepl("year", names)
)

coef_data_alt_spec <- mutate(
  coef_data_alt_spec, 
  names = ifelse(grepl("interc", names, ignore.case = T), "Constant", names),
  names = ifelse(grepl("civil", names, ignore.case = T), "Civil\nwar", names),
  names = ifelse(grepl("neighbor", names, ignore.case = T), "Pct. of\nneighbors", names),
  names = ifelse(grepl("prop", names, ignore.case = T), "Pct. of\nin-degree ties", names),
  names = ifelse(grepl("gle_", names, ignore.case = T), "log of\nGDPC", names),
  names = ifelse(grepl("wdi_", names, ignore.case = T), "Annual\nGDPC growth", names),
  lb = est - 1.96 * se,
  ub = est + 1.96 * se,
  sig = ifelse(ub < 0 | lb > 0, "Statistically\nsignificant","Not\nstatistically significant")
)

ggplot(coef_data_alt_spec, aes(x = est, xmin = lb, xmax = ub, y = reorder(names, -est), col = sig)) + 
  geom_errorbar(width = 0.15, size = 1.25) + 
  theme_bw() + 
  facet_grid(model~dv, space = "free_y", scales = "free_y") + 
  geom_vline(xintercept = 0, size = 1, linetype = "dashed") + 
  labs(x = "Coefficient estimates (95% confidence intervals)", y = "") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 14, color = "black"), 
        axis.text = element_text(size = 12,  color = "black"), 
        axis.title = element_text(size = 14,  color = "black")) + 
  scale_colour_brewer(palette = "Dark2")

ggsave("data_visualization/alt_spec_inst_diff.png", width = 9, height = 10, units = "in")  

load("model_output/unimputed.re.democ.rda")
load("model_output/unimputed.re.autoc.rda")
load("model_output/unimputed.re.rest.part.rda")
load("model_output/unimputed.re.free.media.rda")
load("model_output/unimputed.re.exec.const.rda")

load("model_output/unimputed.glm.democ.rda")
load("model_output/unimputed.glm.autoc.rda")
load("model_output/unimputed.glm.rest.part.rda")
load("model_output/unimputed.glm.free.media.rda")
load("model_output/unimputed.glm.exec.const.rda")

load("model_output/unimputed.lm.democ.rda")
load("model_output/unimputed.lm.autoc.rda")
load("model_output/unimputed.lm.rest.part.rda")
load("model_output/unimputed.lm.free.media.rda")
load("model_output/unimputed.lm.exec.const.rda")

coef_data_no_imputs <- data.frame(
  model = c(
    rep("Random effect logistic", 40),
    rep("Linear probability model", 40),
    rep("Logistic regression", 40)
  ),
  dv = c(
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8)
  ),
  est = c(
    fixef(unimputed.re.democ),
    fixef(unimputed.re.autoc),
    fixef(unimputed.re.exec.const),
    fixef(unimputed.re.free.media),
    fixef(unimputed.re.rest.part),
    
    coef(unimputed.lm.democ),
    coef(unimputed.lm.autoc),
    coef(unimputed.lm.exec.const),
    coef(unimputed.lm.free.media),
    coef(unimputed.lm.rest.part),
    
    coef(unimputed.glm.democ),
    coef(unimputed.glm.autoc),
    coef(unimputed.glm.exec.const),
    coef(unimputed.glm.free.media),
    coef(unimputed.glm.rest.part)
  ),
  names = c(
    names(fixef(unimputed.re.democ)),
    names(fixef(unimputed.re.autoc)),
    names(fixef(unimputed.re.exec.const)),
    names(fixef(unimputed.re.free.media)),
    names(fixef(unimputed.re.rest.part)),
    
    names(coef(unimputed.lm.democ)),
    names(coef(unimputed.lm.autoc)),
    names(coef(unimputed.lm.exec.const)),
    names(coef(unimputed.lm.free.media)),
    names(coef(unimputed.lm.rest.part)),
    
    names(coef(unimputed.glm.democ)),
    names(coef(unimputed.glm.autoc)),
    names(coef(unimputed.glm.exec.const)),
    names(coef(unimputed.glm.free.media)),
    names(coef(unimputed.glm.rest.part))
  ),
  se = c(
    summary(unimputed.re.democ)$coefficients[,2],
    summary(unimputed.re.autoc)$coefficients[,2],
    summary(unimputed.re.exec.const)$coefficients[,2],
    summary(unimputed.re.free.media)$coefficients[,2],
    summary(unimputed.re.rest.part)$coefficients[,2],
    
    summary(unimputed.lm.democ)$coefficients[,2],
    summary(unimputed.lm.autoc)$coefficients[,2],
    summary(unimputed.lm.exec.const)$coefficients[,2],
    summary(unimputed.lm.free.media)$coefficients[,2],
    summary(unimputed.lm.rest.part)$coefficients[,2],
    
    summary(unimputed.glm.democ)$coefficients[,2],
    summary(unimputed.glm.autoc)$coefficients[,2],
    summary(unimputed.glm.exec.const)$coefficients[,2],
    summary(unimputed.glm.free.media)$coefficients[,2],
    summary(unimputed.glm.rest.part)$coefficients[,2]
  )
)

coef_data_no_imputs <- filter(
  coef_data_no_imputs,
  !grepl("year", names)
)

coef_data_no_imputs <- mutate(
  coef_data_no_imputs, 
  names = ifelse(grepl("interc", names, ignore.case = T), "Constant", names),
  names = ifelse(grepl("civil", names, ignore.case = T), "Civil\nwar", names),
  names = ifelse(grepl("neighbor", names, ignore.case = T), "Pct. of\nneighbors", names),
  names = ifelse(grepl("prop", names, ignore.case = T), "Pct. of\nin-degree ties", names),
  names = ifelse(grepl("gle_", names, ignore.case = T), "log of\nGDPC", names),
  names = ifelse(grepl("wdi_", names, ignore.case = T), "Annual\nGDPC growth", names),
  lb = est - 1.96 * se,
  ub = est + 1.96 * se,
  sig = ifelse(ub < 0 | lb > 0, "Statistically\nsignificant","Not\nstatistically significant")
)



ggplot(coef_data_no_imputs, aes(x = est, xmin = lb, xmax = ub, y = reorder(names, -est), col = sig)) + 
  geom_errorbar(width = 0.15, size = 1.25) + 
  theme_bw() + 
  facet_grid(model~dv, space = "free_y", scales = "free_y") + 
  geom_vline(xintercept = 0, size = 1, linetype = "dashed") + 
  labs(x = "Coefficient estimates (95% confidence intervals)", y = "") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 14, color = "black"), 
        axis.text = element_text(size = 12,  color = "black"), 
        axis.title = element_text(size = 14,  color = "black")) + 
  scale_colour_brewer(palette = "Dark2")

ggsave("data_visualization/no_imputs_inst_diff.png", width = 9, height = 11, units = "in")  


load("model_output/gdp.re.democ.rda")
load("model_output/gdp.re.autoc.rda")
load("model_output/gdp.re.rest.part.rda")
load("model_output/gdp.re.free.media.rda")
load("model_output/gdp.re.exec.const.rda")

load("model_output/gdp.glm.democ.rda")
load("model_output/gdp.glm.autoc.rda")
load("model_output/gdp.glm.rest.part.rda")
load("model_output/gdp.glm.free.media.rda")
load("model_output/gdp.glm.exec.const.rda")

load("model_output/gdp.lm.democ.rda")
load("model_output/gdp.lm.autoc.rda")
load("model_output/gdp.lm.rest.part.rda")
load("model_output/gdp.lm.free.media.rda")
load("model_output/gdp.lm.exec.const.rda")


coef_data_gdp <- data.frame(
  model = c(
    rep("Random effect logistic", 40),
    rep("Linear probability model", 40),
    rep("Logistic regression", 40)
  ),
  dv = c(
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8),
    
    rep("Democratic", 8),
    rep("Autocratic", 8),
    rep("Executive\nconstraints", 8),
    rep("Free\nmedia", 8),
    rep("Restricted\npolical\nparticipation", 8)
  ),
  est = c(
    fixef(gdp.re.democ),
    fixef(gdp.re.autoc),
    fixef(gdp.re.exec.const),
    fixef(gdp.re.free.media),
    fixef(gdp.re.rest.part),
    
    coef(gdp.lm.democ),
    coef(gdp.lm.autoc),
    coef(gdp.lm.exec.const),
    coef(gdp.lm.free.media),
    coef(gdp.lm.rest.part),
    
    coef(gdp.glm.democ),
    coef(gdp.glm.autoc),
    coef(gdp.glm.exec.const),
    coef(gdp.glm.free.media),
    coef(gdp.glm.rest.part)
  ),
  names = c(
    names(fixef(gdp.re.democ)),
    names(fixef(gdp.re.autoc)),
    names(fixef(gdp.re.exec.const)),
    names(fixef(gdp.re.free.media)),
    names(fixef(gdp.re.rest.part)),
    
    names(coef(gdp.lm.democ)),
    names(coef(gdp.lm.autoc)),
    names(coef(gdp.lm.exec.const)),
    names(coef(gdp.lm.free.media)),
    names(coef(gdp.lm.rest.part)),
    
    names(coef(gdp.glm.democ)),
    names(coef(gdp.glm.autoc)),
    names(coef(gdp.glm.exec.const)),
    names(coef(gdp.glm.free.media)),
    names(coef(gdp.glm.rest.part))
  ),
  se = c(
    summary(gdp.re.democ)$coefficients[,2],
    summary(gdp.re.autoc)$coefficients[,2],
    summary(gdp.re.exec.const)$coefficients[,2],
    summary(gdp.re.free.media)$coefficients[,2],
    summary(gdp.re.rest.part)$coefficients[,2],
    
    summary(gdp.lm.democ)$coefficients[,2],
    summary(gdp.lm.autoc)$coefficients[,2],
    summary(gdp.lm.exec.const)$coefficients[,2],
    summary(gdp.lm.free.media)$coefficients[,2],
    summary(gdp.lm.rest.part)$coefficients[,2],
    
    summary(gdp.glm.democ)$coefficients[,2],
    summary(gdp.glm.autoc)$coefficients[,2],
    summary(gdp.glm.exec.const)$coefficients[,2],
    summary(gdp.glm.free.media)$coefficients[,2],
    summary(gdp.glm.rest.part)$coefficients[,2]
  )
)

coef_data_gdp <- filter(
  coef_data_gdp,
  !grepl("year", names)
)

coef_data_gdp <- mutate(
  coef_data_gdp, 
  names = ifelse(grepl("interc", names, ignore.case = T), "Constant", names),
  names = ifelse(grepl("civil", names, ignore.case = T), "Civil\nwar", names),
  names = ifelse(grepl("neighbor", names, ignore.case = T), "Pct. of\nneighbors", names),
  names = ifelse(grepl("prop", names, ignore.case = T), "Pct. of\nin-degree ties", names),
  names = ifelse(grepl("gle_", names, ignore.case = T), "log of\nGDPC", names),
  names = ifelse(grepl("wdi_", names, ignore.case = T), "Annual\nGDPC growth", names),
  lb = est - 1.96 * se,
  ub = est + 1.96 * se,
  sig = ifelse(ub < 0 | lb > 0, "Statistically\nsignificant","Not\nstatistically significant")
)



ggplot(coef_data_gdp, aes(x = est, xmin = lb, xmax = ub, y = reorder(names, -est), col = sig)) + 
  geom_errorbar(width = 0.15, size = 1.25) + 
  theme_bw() + 
  facet_grid(model~dv, space = "free_y", scales = "free_y") + 
  geom_vline(xintercept = 0, size = 1, linetype = "dashed") + 
  labs(x = "Coefficient estimates (95% confidence intervals)", y = "") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 14, color = "black"), 
        axis.text = element_text(size = 12,  color = "black"), 
        axis.title = element_text(size = 14,  color = "black")) + 
  scale_colour_brewer(palette = "Dark2")

ggsave("data_visualization/gdp_inst_diff.png", width = 9, height = 11, units = "in")  


load("model_output/autoc_elnet.rda")
load("model_output/democ_elnet.rda")
load("model_output/exec_const_elnet.rda")
load("model_output/free_media_elnet.rda")
load("model_output/rest_part_elnet.rda")
load("model_output/conf_elnet.rda")



autoc_best_tune_coef      <- coef(autoc_elnet$finalModel, autoc_elnet$bestTune$lambda)
democ_best_tune_coef      <- coef(democ_elnet$finalModel, autoc_elnet$bestTune$lambda)
exec_const_best_tune_coef <- coef(exec_const_elnet$finalModel, autoc_elnet$bestTune$lambda)
free_media_best_tune_coef <- coef(free_media_elnet$finalModel, autoc_elnet$bestTune$lambda)
rest_part_best_tune_coef  <- coef(rest_part_elnet$finalModel, autoc_elnet$bestTune$lambda)
conf_best_tune_coef       <- coef(conf_elnet$finalModel, autoc_elnet$bestTune$lambda)


all_coefs <- data.frame(
  names = c(
    names(autoc_best_tune_coef[,1]),
    names(democ_best_tune_coef[,1]),
    names(exec_const_best_tune_coef[,1]),
    names(free_media_best_tune_coef[,1]),
    names(rest_part_best_tune_coef[,1]),
    names(conf_best_tune_coef[,1])
  ),
  est = c(
    autoc_best_tune_coef[,1],
    democ_best_tune_coef[,1],
    exec_const_best_tune_coef[,1],
    free_media_best_tune_coef[,1],
    rest_part_best_tune_coef[,1],
    conf_best_tune_coef[,1]
  ),
  dv = c(
    rep("Autocratic", 8),
    rep("Democratic", 8),
    rep("Constrained executive", 8),
    rep("Free media", 8),
    rep("Restrictions on participation", 8),
    rep("Fatal MIDs", 9)
  )
)

all_coefs <- filter(
  all_coefs, 
  !grepl("Int|year", names)
)

all_coefs <- mutate(
  all_coefs,
  names = as.character(names),
  names = ifelse(grepl("prop", names), "Pct.\nin-degee\nties", names),
  names = ifelse(grepl("neigh", names), "Pct.\nneigh-\nbors", names),
  names = ifelse(grepl("gle_c", names), "log of\nGDPC", names),
  names = ifelse(grepl("wdi_gd", names), "log of\nannual\nGDPC\ngrowth", names),
  names = ifelse(grepl("out_", names), "Out-group\nrelation", names),
  names = ifelse(grepl("civil_", names), "Civil\nwar", names),
  names = ifelse(grepl("democ", names), "Joint\ndemoc-\nracy", names),
  names = ifelse(grepl("autoc", names), "Joint\nautoc-\nracy", names),
  names = ifelse(grepl("major", names), "Joint\nmajor\npowers", names),
  names = ifelse(grepl("cinc", names), "CINC\nratio", names),
  names = ifelse(grepl("contig", names), "Conti-\nguity", names),
  var_to_check = ifelse(grepl("relation|tie", names), "Network\ntopology", "Other\npredictor"),
  dv = factor(
    dv, 
    levels = c(
      "Autocratic",
      "Democratic",
      "Constrained executive",
      "Free media",
      "Restrictions on participation",
      "Fatal MIDs"
    )
  )
)

ggplot(all_coefs, aes(x = reorder(names, -est), y = est, fill = var_to_check)) + 
  geom_col() + 
  facet_wrap(~dv, scales = "free_x", ncol = 2) + 
  theme_bw() + 
  scale_fill_brewer(palette = "Set1") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 14, color = "black"), 
        axis.text = element_text(size = 12,  color = "black"), 
        axis.title = element_text(size = 14,  color = "black")) + 
  labs(x = "", y = "Best tune elastic net estimate")

ggsave("data_visualization/elastic_net_coefs.png", width = 9, height = 11, units = "in")  


## Foreign policy preferences
load("model_output/logit_regression_fp_tau.rda")
load("model_output/logit_regression_fp_kappaba.rda")
load("model_output/logit_regression_fp_pi.rda")
load("model_output/lpm_regression_fp_tau.rda")
load("model_output/lpm_regression_fp_kappaba.rda")
load("model_output/lpm_regression_fp_pi.rda")
load("model_output/logit_fe_regression_fp_tau.rda")
load("model_output/logit_fe_regression_fp_kappaba.rda")
load("model_output/logit_fe_fp_pi.rda")
load("model_output/lpm_fe_fp_tau.rda")
load("model_output/lpm_fe_fp_pi.rda")
load("model_output/lpm_fe_fp_kappaba.rda")


texreg_logit_regression_fp_tau        <- extract.glm(logit_regression_fp_tau)
texreg_logit_regression_fp_kappaba    <- extract.glm(logit_regression_fp_kappaba)
texreg_logit_regression_fp_pi         <- extract.glm(logit_regression_fp_pi)
texreg_lpm_regression_fp_tau          <- extract.glm(lpm_regression_fp_tau)
texreg_lpm_regression_fp_kappaba      <- extract.glm(lpm_regression_fp_kappaba)
texreg_lpm_regression_fp_pi           <- extract.glm(lpm_regression_fp_pi)
texreg_logit_fe_regression_fp_tau     <- extract.fixest(logit_fe_regression_fp_tau)
texreg_logit_fe_regression_fp_kappaba <- extract.fixest(logit_fe_regression_fp_kappaba)
texreg_logit_fe_regression_fp_pi      <- extract.fixest(logit_fe_fp_pi)
texreg_lpm_fe_fp_tau                  <- extract.fixest(lpm_fe_fp_tau)
texreg_lpm_fe_fp_pi                   <- extract.fixest(lpm_fe_fp_pi)
texreg_lpm_fe_fp_kappaba              <- extract(lpm_fe_fp_kappaba)


texreg(list(texreg_logit_regression_fp_tau,
            texreg_logit_regression_fp_kappaba,
            texreg_logit_regression_fp_pi,
            texreg_lpm_regression_fp_tau,
            texreg_lpm_regression_fp_kappaba,
            texreg_lpm_regression_fp_pi,
            texreg_logit_fe_regression_fp_tau,
            texreg_logit_fe_regression_fp_kappaba,
            texreg_lpm_fe_fp_pi,
            texreg_lpm_fe_fp_tau,
            texreg_lpm_fe_fp_kappaba,
            texreg_lpm_fe_fp_pi), stars = 0.05, digits = 2)


si_table_8 <- texreg(list(texreg_logit_regression_fp_tau,
                          texreg_logit_regression_fp_kappaba,
                          texreg_logit_regression_fp_pi,
                          texreg_lpm_regression_fp_tau,
                          texreg_lpm_regression_fp_kappaba,
                          texreg_lpm_regression_fp_pi,
                          texreg_logit_fe_regression_fp_tau,
                          texreg_logit_fe_regression_fp_kappaba,
                          texreg_logit_fe_regression_fp_pi,
                          texreg_lpm_fe_fp_tau,
                          texreg_lpm_fe_fp_kappaba,
                          texreg_lpm_fe_fp_pi), stars = 0.05, digits = 2)

write.table(si_table_8, file = "data_visualization/si_table_8.txt")
